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Abstract 


Our  goal  is  the  design  of  fast  parallel  algorithms  statistical  inference  for  scalar 
and  wave  propagation  problems.  We  have  looked  at  source  inversion  and 
inverse  medium  problem  problems.  We  use  a  Bayesian  approach  in  which  the 
regularization  appears  as  prior  information  and  the  data  mismatch  appears  as  a 
likelihood  information,  given  known  noise  probability  density  functions. 

A  key  component  of  all  of  our  algorithms  is  the  approximation  of  the  Hessian 
operator.  Key  components  of  our  work  are  rank-revealing  factorizations,  fast 
extraction  of  the  diagonal  of  the  inverse,  adaptivity,  and  integration  of  all  of  these 
components  within  a  particle  filter  methodology. 


In  addition,  our  implementations  are  being  designed  to  scale  on  manycore  and 
heterogeneous  parallel  architectures. 


Figure  1:  To  demonstrate  the  effectiveness  of  our  algorithms  we  have  solved  a  multipoint  illumination  inverse 
scattering  problem.  Here  the  scatterer  is  an  aircraft  (Lockheed  SR-71 ).  The  first  figure  from  left  depicts  the 
discretization  mesh;  the  second  figure  depicts  the  illumination  due  to  the  incoming  field;  the  third  figure  depicts  the 
scattered  field  (data)  on  a  detector  plane  placed  at  a  distance  below  the  aircraft;  the  fourth  image  depicts  a 
reconstruction  of  the  scatterer  using  the  algorithms  developed  in  our  group.  These  simulations  were  done  using  the 
Born  approximation  for  the  scalar  Helmholtz  equation  using  a  boundary  integral  equation  formulation  accelerated 
using  our  kernel-independent  Fast  Multipole  Method.  In  the  inversion  problem,  we  seek  to  combine  multiple  frequency 
and  multiple  direction  illumination  data  (rows  two  and  four)  to  reconstruct  the  scatterer  solving  the  inverse  medium 
problem.  We  give  more  details  on  the  reconstruction  in  figures  2  and  3. 


Introduction 


Our  goal  is  to  design  scalable  algorithms  for  inverse  medium  and  inverse  source 
problems  in  wave  propagation  for  both  time-harmonic  and  time-domain 
formulations. 

In  inverse  scattering  problems,  we  seek  to  reconstruct  infinite-dimensional  fields 
by  combining  models  of  the  underlying  physics  with  observations;  typically  the 
underlying  input-output  (parameter-to-observable)  operators  are  compact  and 
result  in  ill-posed  inverse  problems.  Severe  computational  costs  of  existing 
methods  (e.g.,  associated  with  multiple  forward  and  adjoint  problem  solutions, 
large  dense  operators,  and  strong  nonlinearities  even  when  the  forward  problem 
is  linear)  limit  problem  size  and  model  complexity.  In  addition,  designing  effective 
solution  methods  that  scale  well  requires  application  of  problem-specific. 

In  such  problems,  one  is  given  a  noisy  measurement  of  the  scattered  wave  at 
certain  receiver  locations  and  seeks  to  reconstruct  a  medium  inhomogeneity  or 
an  unknown  source.  The  forward  problem  can  be  written  as  d  =  K  z,  where  z  is 
the  inversion  parameter,  K  is  the  input-ouput  operator,  and  d  are  the 
measurements  in  the  receivers.  The  deterministic  inverse  problem  is  to 
reconstruct  z  given  d,  that  is  invert  K.  The  probabilistic  inverse  problem  is  to  get 
the  probability  density  function  of  z,  given  some  prior  knowledge  and  known 
noise  characteristics. 

The  operator  K  plays  a  key  role  in  both  deterministic  and  probabilistic  inversions. 
The  Hessian  H  is  defined  as  H=KTK. 

Applying  K  to  a  vector  involves  the  solution  of  wave  propagation  problems  and 
applying  KT  involves  the  adjoint  operator.  In  fact,  for  problems  that  involve 
multiple  illuminations  and  multiple  frequencies  then  a  matrix-vector  multiplication 
with  the  Hessian  involves  numerous  forward  and  adjoint  wave  propagations. 
Finally,  the  Hessian  is  typically  a  compact  operator  and  thus,  is  not  invertible  in  a 
stable  manner. 

Status/Progress 

In  this  report  we  discuss  progress  on  several  related  fronts  together  aimed  at 
developing  and  applying  fast  scalable  methods  for  quantification  of  uncertainty  in 
the  solution  of  inverse  wave  propagation  problems: 

•  FalMS:  It  is  a  fast  approximation  scheme  for  time-harmonic  problems  with 
multiple  illuminations.  But  it  can  be  integrated  with  the  time-domain 
solvers  developed  by  Prof.  Ghattas’s  group.  In  fact  the  solver  can  be 
used  for  any  problem  in  which  two  properties  are  necessary,  the  existence 
of  a  fast  method  to  evaluate  the  Green’s  function  of  the  underlying 
differential  operator,  and  smoothness  of  this  Green’s  function,  which  is 
required  in  order  to  compress  long-range  interactions. 


•  Adaptivity  for  inverse  medium  problems:  We  are  developing  algorithms 
specific  to  the  inverse  medium  problem.  We  are  exploring  uncertainty- 
based  adaptivity  and  jump-localization  techniques.  The  hardest  problem  in 
adaptivity  is  the  interplay  between  the  regularization  and  the  discretization. 
Indeed,  discretization  corresponds  to  some  form  of  regularization.  We 
have  developed  a  scheme  that  shows  promising  results  for  simple 
examples. 

FalMS. 

Our  recent  contribution  is  the  development  of  FalMS,  and  algorithm  for  the 
inverse  medium  problem  for  the  scalar  wave  equation.  We  consider  the  inverse 
medium  problem  for  the  time-harmonic  wave  equation  with  broadband  and  multi¬ 
point  illumination  in  the  low  frequency  regime.  Such  a  problem  finds  many 
applications  in  geosciences  (e.g.  ground  penetrating  radar),  non-destructive 
evaluation  (acoustics),  and  medicine  (optical  tomography).  We  use  an  integral- 
equation  (Lippmann-Schwinger)  formulation,  which  we  discretize  using  a 
quadrature  method.  We  consider  only  small  perturbations  (Born  approximation). 
To  solve  this  inverse  problem,  we  use  a  least-squares  formulation.  We  present  a 
new  fast  algorithm  for  the  efficient  solution  of  this  particular  least-squares 
problem.  Let  us  emphasize  that 


Figure  2:  In  this  figure,  we  demonstrate  the  problem  of  a  multipoint  illumination  inverse  scattering  problem.  In  the  first 
row  (from  the  top),  we  depict  the  plane-wave  illumination  of  the  scatterer  for  different  frequencies  (which  correspond  to 
a  scatterer  size  equal  to  1,4,16, and  64  wavelengths  respectively).  In  the  second  row,  we  show  the  corresponding 
amplitude  of  the  scattered  field  due  to  the  illumination.  The  observations  are  made  on  a  plane  at  a  distance  below  the 
scatterer  and  detector  resolution  is  sufficient  to  resolve  the  minimum  wavelength.  In  the  third  row,  we  show  multiple 


illuminations  using  the  same  wavelength  but  different  illumination  directions;  in  the  fourth  row  we  show  the 
corresponding  scattered  field  on  the  detector  plane.  The  results  of  the  inversion  are  given  in  Figure  3. 


If  Nf  is  the  number  of  excitation  frequencies,  Ns  the  number  of  different  source 
locations  for  the  point  illuminations,  Nd  the  number  of  detectors,  and  N  the 
parametrization  for  the  scatterer,  a  dense  singular  value  decomposition  for  the 
overall  input-output  map  will  have  [min(Ns  Nf  N,  N)]A{2}  X  max(Ns  Nf  Nd,  N) 
cost.  We  have  developed  a  fast  SVD-based  preconditioner  that  brings  the  cost 
down  to  0(  (N  +  Ns  +  Nf+  Nd)  log  N  )  thus,  providing  orders  of  magnitude 
improvements  over  a  black-box  dense  SVD  and  an  unpreconditioned  linear 
iterative  solver. 


Figure  3:  In  this  figure,  we  depict  reconstructions  for  different  illumination  and  maximum  frequency  regimes  using  our 
FAIMS  inverse  medium  solver.  We  use  data  created  by  the  simulations  and  the  scatterer  (SR-71 )  described  in  Figure  2. 
The  scatterer  is  modeled  as  an  unknown  medium  perturbation  (point  scatterers).  In  the  top  row  from  left,  the  first 
reconstruction  was  done  using  a  single  sub-wavelength  illumination,  the  second  using  an  incoming  wave  whose 
frequency  corresponds  to  the  true  scatterer  being  one  wavelength  long  and  8  illuminations,  the  third  using  a  four- 
wavelength  frequency  and  8  illuminations,  and  the  last  one  using  a  four-wavelength  frequency  and  32  illuminations.  In 
the  second  row  we  depict  different  views  (top, side, front, and  angled)  of  a  high  resolution  reconstruction  using  four 
different  illumination  frequencies  (the  maximum  illumination  frequency  corresponds  to  eight  wavelengths)  and  128 
illumination  directions  (for  each  frequency).  As  we  increase  the  illumination  frequency,  the  rank  of  the  Hessian 
increases  and  the  reconstructions  become  computationally  demanding.  Nevertheless,  we  are  able  to  reconstruct  the 
target  scatterer  quite  accurately. 


This  work  has  been  completed  and  a  paper  has  been  submitted  to  SISC.  Our 
next  step  is  to  extend  it  to  nonlinear  inverse  medium  problems  and  to  introduce 
adaptivity.  Also,  we  are  developing  a  parallel  library  that  implements  FalMS.  In 
the  application  domain,  we  are  extending  FalMS  to  the  Maxwell  and  Navier 
equations  (time-harmonic)  and  we  are  integrating  with  the  massively  parallel  Fast 
Multipole  Methods  developed  in  our  group. 

Adaptivity. 

We  are  investigating  adaptive  methods  for  the  inverse  medium  for  the  Helmoltz 
equation.  The  medium  is  represented  using  an  octree  data-structure  and 
piecewise  Chebyshev  polynomials.  We  are  investigating  two  approaches,  one  is 
based  on  the  Hessian  information  and  the  second  one  is  based  on  Tadmor 
concentration  kernels. 

In  the  first  approach,  the  diagonal  of  the  inverse  of  the  Hessian  gives  us  the 
diagonal  of  the  covariance  matrix,  which  in  turn  indicates  the  uncertainty 
(variance)  of  our  estimate.  If  this  uncertainty  is  high,  we  do  not  refine  otherwise 


we  do.  The  diagonal  of  the  inverse  we  will  be  built  by  a  randomization  method 
that  only  requires  matrix-vector  multiplications, 
detectors 


5  A 

Figure  4:  LEFT:  Settings  of  the  simplified  configuration  considered  in  this  example.  The  detectors  (blue  circles)  are 
located  on  a  plane  (Nd  =  100).  The  sources  (red  circles)  are  located  on  a  sphere  (Ns  =  12).  We  assume  a  line  singularity  of 
the  scatterer,  with  strength  given  by  the  inset  at  the  lower  left.  We  use  Chebyshev  polynomials  for  the  discretization  of 
the  forward  problem,  along  with  a  direct  evaluation  of  the  The  discretization  is  ID  (along  the  dashed  line).  The  number 
of  frequencies  is  set  to  Nf  =3. The  size  of  the  box  is  5X.  RIGHT:  We  report  results  for  our  adaptive  scheme  in  the  absence  of 
noise. 

In  the  second  approach,  is  computationally  cheaper  as  it  is  entirely  local.  It  uses 
edge  detection  using  the  spectral  coefficients  of  a  discontinuous  function.  This 
estimator  will  be  based  on  the  work  of  Anne  Gelb  and  Eitan  Tadmor,  Acta 
Numerica  1999.  The  basic  idea,  is  that  upon  inversion  the  Chebyshev 
coefficients  of  every  octree  leaf  are  examined  with  the  help  of  a  concentration 
kernel,  and  if  the  resulting  jump  is  sufficiently  large  (this  is  an  application- 
dependent  parameter),  then  we  proceed  with  refinement.  This  process  will  be 
combined  with  multiscale  sampling  techniques  for  Bayesian  inversion  schemes. 

Eventually,  the  two  adaptivity  approaches  will  be  combined  as  the  uncertainty- 
based  one  does  not  use  information  from  the  reconstructed  field  and  the  jump- 
detection  one  does  not  address  the  ill-posedness  of  the  inverse  problem.  The 
work  has  not  been  completed  yet,  but  we  continue  our  work  with  funding  by  other 
resources.  Preliminary  results  are  reported  below. 


Acknowledgment/Disclaimer 

This  work  was  sponsored  (in  part)  by  the  Air  Force  Office  of  Scientific  Research, 
USAF,  under  grant/contract  number  FA9550-09-1-0608.  The  views  and 


conclusions  contained  herein  are  those  of  the  authors  and  should  not  be 
interpreted  as  necessarily  representing  the  official  policies  or  endorsements, 
either  expressed  or  implied,  of  the  Air  Force  Office  of  Scientific  Research  or  the 
U.S.  Government. 

Personnel  Supported 

Stephanie  Chaillat 
George  Biros 

Publications 

•S.  Chaillat  and  G.  Biros,  An  adaptive  algorithm  for  the  nonlinear  inverse 
medium  problem  for  the  Helmholtz  equation,  in  preparation 
•S.  Chaillat  and  G.  Biros,  FalMS:  A  fast  algorithms  for  the  inverse  medium 
problem  with  multiple  frequencies  and  multiple  sources  for  the  scalar  Helmholtz 
equation.  Journal  of  Computational  Physics,  231  (20),  2012 
•S.S.  Adavani  and  G. Biros,  Fast  algorithms  for  inverse  problems  with  parabolic 
PDE  constraint,  SIAM  Journal  on  Imaging  Sciences,  SIAM  Journal  on  Imaging 
Sciences,  to  appear 

•S.S.  Adavani  and  G.  Biros,  Fast  algorithms  for  source  identification  problems 
with  elliptic  PDE  constraints,  SIAM  Journal  on  Imaging  Sciences,  3  (4),  2010 
•E.G.  Mendizabal-Ruiz,  G.  Biros  and  I.  Kakadiaris,  An  inverse  scattering 
algorithm  for  the  segmentation  of  the  luminal  border  on  intravascular 
ultrasound  data,  Proceedings  of  MICCAI  2009 
•  I.  Lashuk,  G.  Biros  et  al,  A  massively  parallel  adaptive  fast-multipole  method  on 
heterogeneous  architectures,  Proceedings  of  SC  2009 


During  Duration  of  Grant 

Research  Associate,  Georgia  Tech 
Associate  Professor,  Georgia  Tech 
Professor,  The  University  of  Texas  at  Austin 


Honors  and  Awards  Received 

•The  paper  I.  Lashuk,  A.  Rahimian,  G.  Biros  et  al.  “Petascale  direct  numerical 
simulation  of  blood  flow  on  200K  cores  and  heterogeneous  architectures,” 
received  the  2010  Gordon  Bell  Prize. 

•The  paper  I.  Lashuk,  G.  Biros  et  al,  A  massively  parallel  adaptive  fast-multipole 
method  on  heterogeneous  architectures,  Proceedings  of  SC  2009,  was  a  Best 
Technical  Paper  finalist. 

•Biros  gave  a  plenary  talk  on  fast  elliptic  solvers  at  the  2010  Parallel  Matrix 
Algorithms  and  Applications  Conferences,  in  Basel,  Switzerland,  July  2010 

•  Biros  gave  a  plenary  talk  on  inverse  problems  at  the  Applied  Inverse  Problems, 
Conference  in  Vienna,  July  2009 


AFRL  Point  of  Contact:  None 
Transitions:  None 

New  Discoveries:  Fast  algorithms  for  inverse  medium  solvers 


